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INTRODUCTION 


This  is  the  first  semiannual  technical  report  under  Contract 
F30602-74-C-01 3CT enti tied  Atmospheric  Effects  on  Space  Object  Imagery. 

The  report  covers  efforts  during  the  period  January  15,  1974  to  June  1 
1974.  One  object  of  the  program  is  to  study  techniques  for  the  computer 
generation  of  random  wavefronts  simulating  atmospherically  degraded 
light  beams.  This  material  is  intended  to  be  useful  in  the  computer 
simulation  of  restoration  and  predetection  compensation  of  atmospherically 
degraded  optical  images,  especially  images  of  high  flying  objects  such 
as  satellites.  A second  objective  is  to  study  angle  of  arrival  fluctu- 
ations and  high  altitude  temperature  spatial  spectra.  This  work  will  be 
considered  in  a later  report. 

In  the  design  and  construction  of  equipment  to  overcome  the  effects 
of  atmospheric  turbulence,  it  is  often  a difficult  procedure  to  cneck  out 
the  equipment  in  the  environment  for  which  it  is  intended.  Large 
telescopes  are  often  busy  and  ground-level  propagation  ranges  are  not 
always  accessible.  For  these  reasons  it  is  desirable  to  simulate  the 
operation  of  such  devices  digitally.  In  order  to  perform  digital  simu- 
lation, one  must  have  as  exact  a representation  of  the  turbulent  atmosphere 
as  possible.  This  report  concerns  efforts  to  provide  efficient  computer 
simulation  of  atmospherically  degraded  light  beams  as  a first  step  in 
simulation  of  a complete  system. 

Some  information  is  available  for  simulating  atmospherically  de- 
graded wavefronts.  It  is  believed  that  the  phase  variations  follow 
approximately  Gaussian  statistics  as  do  phase  difference  fluctuations 
[Clifford,  1971]  [Bertolotti , 1968]  and  log-amplitude  variations 
[Ochs,  1969],  [Tatarski,  1971,  p.  292],  It  is  further  expected  that  the 
phase  structure  function  will  follow  a +5/3  law  with  separation  in  the 
inertial  subrange  [Bouricius,  1970],  Further,  predicted  expressions  for 
phase  structure  and  correlation  functions  are  known  [Tatarski,  1971]. 

These  items  are  sufficient  to  allow  the  generation  of  a set  of  randomly 
degraded  wave  fronts. 

In  the  present  contract  the  objectives  are  to  simulate  random 
optical  fields  using  two  approaches,  a random  matrix  approach  and  the 
Karhunen-Loeve  series  approach,  and  to  examine  the  statistics  of 
atmospherically  degraded  images  using  recently  obtained  stellar  image 
data.  The  random  matrix  approach  has  been  used  by  others  [McGlamery, 

1974],  [Bradley,  1974]  to  generate  random  phase  fronts,  but  the 
processes  of  generating  the  associated  log-amplitude  fluctuations  and 
of  simulating  the  fields  for  a light  beam  that  has  propagated  downward 
along  a vertical  path  have  not  been  hitherto  considered.  Neither  has 
the  application  of  the  Karhunen-Loeve  series  to  this  problem  been 
considered.  During  the  past  quarter  our  efforts  have  concentrated  on 
setting  up  and  examining  matrices  of  random  numbers  to  assure  that 
they  satisfy  the  statistical  criteria.  Work  on  the  Karhunen-Loeve 
approach  is  just  starting.  The  work  on  the  examination  of  stellar 
images  has  not  yet  ctarted. 
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In  the  present  report  the  work  on  the  generation  of  random  -‘base 
matrices  will  be  described  first,  followed  by  a derivation  of  the 
scheme  for  relating  phase  and^amplitude  fluctuations.  Finally  a b'ief 
discussion  of  the  Karhunen-Loeve  approach  will  follow. 

The  procedure  for  generating  matrices  of  random  numbers  with  the 
desired  statistics  is  quite  straightforward.  The  fmt  step  is  to 
generate  an  ensemble  of  arrays  of  random  numbers  with  normal  distri- 
butions and  uniform  average  spectra.  The  second  step  is  to  compute 
the  spectrum  of  each  individual  array,  tailor  these  spectra  to  the 
desired  average  functional  dependence,  and  transform  back  to  obtain 
random  arrays  with  the  desired  probability  density  and  covariance. 

Some  of  the  details  of  the  process  we  used  will  now  be  presented. 
The  first  step  was  to  obtain  a set  of  arrays  of  uncorrelated  random 
numbers.  These  random  number  arrays  consisted  of  matrices  of  64  x 64 
elements  which  were  (theoretically)  samples  from  a Gaussian  probability 
distribution.  These  elements  were  produced  by  a machine-specific  random 
number  generator. 

To  insure  that  there  be  no  correlation  effects  between  the  elements 
of  the  array  a random  loading  algorithm  was  employed.  In  this  algorithm 
two  independent  random  number  generators  were  used  to  generate  numbers 
between  one  and  sixty-four.  These  random  number  tuples  were  used  to 
choose  a point  in  the  array  which  was  then  filled  by  a third  generator 
with  Gaussian  random  numbers  between  -3  and  +3.  This  process  was 
continued  until  the  array  was  something  more  than  half-filled.  The  re- 
maining cells  were  then  filled  sequentially. 

The  random  loading  scheme  just  described  obviously  has  some  degr  : 
of  redundancy  due  to  the  possibility  of  selecting  the  same  point  in  the 
array  more  than  once,  but  in  the  interest  of  generating  arrays  in  which 
there  is  no  element  interaction  (correlation),  this  redundancy  was  felt 
worthwhile. 

The  arrays  were  then  tested  for  the  desired  statistical  properties. 
This  was  accomplished  by  grouping  the  arrays  by  fives  and  doing  an 
analysis  of  variance  with  two-way  classification  on  each  group.  This 
test  assumes  a model  of  the  form  shown  in  Eq.  (1). 

O)  xij  = w + «i  + 

where  x-ji  is  the  generic  array  element;  v is  the  mean;  ai,  6j,  and  aij 
represent  respectively  row  effects,  column  effects,  and  row-column 
interactions  or  correlations.  Statistics  were  generated  from  the  arrays 
enabling  us  to  perform  F-tests  for  significance  of  these  effects.  No 
such  significant  effects  were  found. 
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The  next  step  was  to  obtain  the  spatial  spectrum  of  the  Individual 
arrays  and  to  adjust  the  power  spectra  to  have  the  desired  functional 
dependence.  This  was' accomplished  using  the  Fast  Fourier  Transform  (FFT) 
technique.  After  the  spectral  components  were  found,  they  were  put  in 
phasor  form  giving  the  magnitude  and  phase  of  each  component.  The  mag- 
nitude was  then  multiplied  by  the  value  of  the  desired  phase  spatial 
spectrum  for  that  value  spectral  component.  The  phase  spatial  spectrum 
chosen  was  based  on  the  Von  Karman  index  spectrum  with  a ten  meter  outer 
scale.  This  was  available  from  work  performed  previously  but  was 
normalized  to  coherence  length,  r0,  as  discussed  In  Appendix  B.  An 
inner  scale  much  smaller  than  the  aperture  grid  size  was  assumed.  The 
complex  phase  of  each  spectral  component  was  saved  for  use  later  In 
generating  the  log-amplitude.  The  modified  spectra  were  then  transformed 
back  to  the  spatial  domain.  The  result  was  the  desired  set  of  random 
arrays  with  the  proper  phase  correlation  function  and  normal  probability 
distribution  function. 

After  the  phase  fronts  were  generated  and  checked  statistically, 
the  associated  log-amplitude  fluctuations  were  considered.  The 
derivation  of  the  procedure  for  relating  the  log- amplitude  to  the  phase 
will  now  be  presented.  The  discussion  uses  the  well-known  complex  log- 
amplitude  representation.  We  assume  that  the  (random)  phase  of  the 
light  wave  entering  a receiving  telescope  is  known,  and  that  the  region 
inside  the  telescope  has  no  index  fluctuations  so  that  the  free  space 
wave  equation  applies,  e*j(r)  = 0.  The  complex  log-amplitude  4>(r)  1s» 
following  standard  procedure,  divided  into  an  unperturbed  plane  wave, 

-jkz,  with  constant  amplitude,  +Aq»  and  a correction,  f|»  due  to  the  random 
fluctuation  as  in  Eqs.  (2).  The  correction  follows  the  differential 
equation  in  Eq.  (3)  [Tatarski,  1971,  p.  223]. 

(2a)  E(r)  = e*^ 


(2b)  i|i  ■ - jkz  + 
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Using  Eq.  (2c)  we  can  separate  Eq.  (3)  into  real  and  imaginary  parts 
giving  two  equations  relating  ^(r)  and  f|(r).  They  are 
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(4b)  = -2  k If- 

ax*  ay*  92  * 

The  two  coupled  differential  equations  relating  the  log-ampl1*;>de 
and  phase  are  easily  solved.  Solutions  are  given  In  Eqs.  (5). 

(5.)  ♦ ■ S «(*«  ♦ ♦ (|2? )*  ♦ <^)2)k  * e„) 


(5b)  . - n *„  bin(^  * ♦ fep )*  . (^n.)  j If  * ,„] 

They  are  In  Fourier  series  form  to  facilitate  matching  to  the  computer 
generated  random  phase  and  log-anplltude  arrays.  We  note  that  the 
relationship  between  the  phase  and  log-amplitude  expressions  Is 
especially  simple;  the  only  difference  Is  that  cosines  are  replaced  by 
sines. 

Equations  (5)  are  more  than  sufficient  for  fitting  the  computer- 
generated solutions.  The  summations  cover  the  range  -«xm,n<«  while  only 
the  range  <m,n<«  Is  required  to  fit  a given  random  matrix;  l.e.,  two 
Independent  solutions  are  possible.  Two  such  solutions  are  giver  In 
Eqs.  (6)  and  (7).  They  differ  by  choice  of  the  relations  between  the 
Independent  constants. 


(6a) 

0 = - 0 
mn  -m,-n 

(6b) 

* - + a 

ym,n  -m,-n 

(7a) 

em,n  3 " 9-m,-n 

(7b) 

A m - < 

Tm,n  -m,-n 

The  solutions  represented  by  Eqs.  (6)  and  (7)  are  chosen  for  a specific 
reason.  That  represented  by  Eqs.  (5)  and  (6)  has  finite  phase  and  zero 
log-amplitude  In  the  plane  z=0.  The  solution  represented  by  Eqs.  (5) 
and  (7)  have  zero  phase  and  finite  log-amplitude  at  z*0.  The  sum  of 
the  two  solutions  can  be  used  to  represent  any  Independently  generated 
combination  of  log-amplitude  and  phase. 
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The  procedure  thus  suggested  is  to  generate  two  statistically 
independent  random  fields  for  log-amplitude  and  phase  respectively.  This 
procedure  is  especially  appropriate  for  a plane  wave  propagated  hori- 
zontally through  the  atmospnere  a sufficiently  long  distance,  L»2ira2/x; 
(a  is  the  receiver  aperture  size).  For  such  ranges  the  cross  covariance 
Bju(p)  is  negligibly  small  [Tatarski,  1971,  Eqs.  (36)-(38),  (46)].  How- 
ever, the  assumption  of  zero  cross  covariance  for  vertical  paths  where 
the  effective  range  is  limited  is  not  so  obvious. 

A second  technique  being  investigated  for  generating  the  random 
phase  fronts  involves  the  use  of  the  Karhunen-Loeve  series  [Davenport, 
1958].  Here  the  correlation  function  B(r,r')  is  automatically  built 
into  the  functions  because  it  is  the  kernal  of  the  integral  equation 

(6)  | B(r,r')  fn(r)  d r = *n  fn(r  ) 

used  to  generate  {f«},  the  set  of  orthonormal  functions.  There  is 
further  only  a one-dimensional  set  of  random  numbers  to  consider,  namely 
the  random  function  coefficients  in  the  series  for  the  phase  front.  One 
great  advantage  in  the  use  of  the  K-L  series  is  the  combined  statistical 
as  well  as  functional  orthogonality  property  of  the  functions. 

In  the  Karhunen-Loeve  procedure  the  phase  covariance  derived 
from  a Von  Karmann  index  spatial  spectrum  was  used.  An  8 x 8 element 
array  square  aperture  gave  sufficient  precision  for  the  lower  order  poly- 
nomials. The  procedure  is  quite  similar  to  one  used  previously 
[Moreland,  1969]. 

The  procedure  for  generating  the  set  of  random  phase  fronts  is 
then  to  generate  the  series  of  random  coefficients,  an,  for  each  member 
of  the  set  using  a random  number  generator  and  sum  at  every  point  in 
the  array.  Once  the  phase  fronts  have  been  generated,  the  random  log- 
amplitude  functions  are  generated  as  indicated  previously. 

This  procedure  has  the  advantage  that  only  a one  dimensional  set 
of  random  numbers  needs  to  be  generated  per  array,  but  has  the  dis- 
advantage that  the  series  must  be  summed  for  every  point  in  every  member 
of  every  set. 

To  summarize,  two  methods  of  generating  random  phase  and  log- 
amplitude  fields  are  being  investigated.  The  first  method  generates 
a set  of  random  arrays  and  tailors  them  so  as  to  have  the  requisite 
probability  distribution  and  covariance.  The  second  method  generates 
a set  of  eigenfunctions  having  the  correct  covariance.  The  series 
expansions  for  the  phase  and  log-amplitude  are  then  summed  using  random 
coefficients  to  provide  one  member  of  the  set  of  wavefronts.  We  are 
also  studying  the  relationship  between  the  random  arrays  and  atmos- 
pherically degraded  log-amplitude  and  phase  fronts  to  assure  that  the 
computer  generated  arrays  represent  as  accurately  as  possible  the  real 
atmosphere. 
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APPENDIX  A 


In  this  appendix  the  random  number  generation  scheme  mentioned 
in  the  text  is  described.  The  random  number  generators  are  based  on 
the  power  residue  method  (IBM  6C20-8011-0).  This  method  consists  of  the 
following  steps: 

1,  Choose  for  a starting  value  any  odd  integer  yQ. 

2,  Choose  a constant  multiplier  of  the  form  x = 8t±3  where  t is 
an  integer  such  that  x is  close  to  2b/2  where  the  machine 
word  length  is  b bits. 

3,  Compute  the  product  xu0.  This  product  is  2b  bits  long.  The 
higher  order  b bits  are  discarded  and  the  lower  order  b bits 
become  the  next  number  u-j . 

In  this  manner  the  random  numbers  are  generated  according  to  the  formula 

Vi  s xv 

It  is  understood  in  the  above  discussion  that  the  radix  poim.  in 
the  u-'s  is  at  the  extreme  right  of  the  word  so  that  the  nur.^ers  are 
integers.  After  generation  however  the  radix  point  is  moved  to  the 
extreme  left  of  the  word  (by  dividing  by  the  maximum  representable  integer) 
to  produce  a number  between  zero  and  unity.  At  this  point  the  random 
nuiribers  can  be  scaled  to  produce  the  required  range  of  values,  e.g.,  1-64. 


APPENDIX  B 


In  this  appendix  we  present  the  derivation  of  the  expression  for 
the  spatial  spectrum  of  the  phase  structure  function.  The  situation 
considered  has  a plane  wave  propagating  vertically  downward  through  the 
earth 's  boundary  layer  where  the  turbulence  structure  constant  is  a 
function  of  altitude.  The  outer  scale  is  also  a function  of  altitude 
but  this  will  be  neglected  because  it  has  little  effect  on  the  resulting 
expressions.  The  object  is  to  present  an  expression  for  the  phase 
structure  function  spectrum  normalized  in  a simple  fashion;  in  this  case 
normalized  to  the  appropriate  value  for  the  coherence  length,  r . 

The  basic  expressions  for  obtaining  the  phase  struction  function 
spectrum  have  been  given  [Tatarski,  1971].  Combining  Eqs.  (8)  and  (9) 
on  p.  243  and  Eq.  (15a)  oji  p.  230  and  substituting  for  the  relative 
permittivity  spectrum  $e(ic)  in  terms  of  the  refractive  index  spectrum, 
as  in  Eq.  (Bl), 

(Bl)  *c(k)  = 4$n(x)  - 4 cjj(n)  *noM> 

we  have  for  the  plane  wave  phase  structure  function  spectrum 

(B2a)  Fs(k)  = *k2  *no(K,0)  £ (l  + C2(K)  dn  . 

The  corresponding  expressions  for  the  phase  structure  function  and 
covariance  are  respectively 

(B2b)  D$(p)  = 4tt  f (1  - Jq(kp))  F$(k)  <d< 


(B2c)  B$(p)  ■ 2ir  I JQ  (kp)  F$(k)  <dK  . 


In  Eqs.  (1)  and  (2)  n is  the  longitudinal  variable  running  from  trans- 
mitter to  receiver  and  k is  the  magnitude  of  the  transverse  components 

of  spatial  frequency,  * =J  + *(■  . (k  and  n are  in  the  longitudinal 
direction.)  v x y z 

For  spatial  frequencies  corresponding  to  fluctuations  in  the  size 
range  5 mm  to  5 m,  for  ranges  in  excess  of  1 km,  and  for  visible  light 
we  can  replace  the  cosine  term  in  Eq.  (B2a)  by  its  average,  zero.  We 
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further  use  the  height  variation  of  the  index  structure  parameter,  C2(h) 
for  well  developed  turbulence  in  unstable  air  [Wyngaard,  1971] 

o o u **4/3 
(S3)  Cp(h)  - C2(h0)  (J-) 

0 

Thus,  taking  the  source  at  a very  large  height  and  receiver  near  ground 
level  we  have 

2 (L  ? / H +L-n  \‘4/3 

(B4)  FS<K>  * *no^<,0)  Jo  Cn  (Ho)(— ) dn 


o 

where  H0  is  the  height  of  the  receiver  and  Cn(H0)  is  the  structure 

parameter  value  at  the  receiver  height,  (n^L) - Performing  the  n inte- 
gration and  substituting  for  $no(K»°)  using  the  Von  Karmann  spectrum 
we  have 


(B5) 


Fc(k)  = ir<2*3H^(Hj[l-(l  + ±-) 


-1/3' 


o n'  o 


x.033 


r 2 r11/6 

(L»ZI, 

0 


Taking  the  source  much  higher  than  the  receiver,  L + H0  » H0  and 
rearranging,  we  have 


(B6) 


Fs(k)  » .033Trk2  (3Ho)C2(Ho) 


1.077/ 

~r~> 


-11/6 


+ K 


where  L0  is  the  outer  scale  and  the  factor,  1.077,  is  included  to  make 
the  asymptotic  value  of  the  index  structure  for  large  separation, 

Dn(")  - This  allows  convenient  determination  of  U from  a 

log-log  plot  of  Dn  versus  p:  the  break  point  occurs  where  the  separation 
equals  L0. 

Equation  (6)  is  the  desired  spectrum.  However  it  will  be  more 
convenient  to  express  it  using  the  coherence  length,  r0,  for  normalization. 
For  this  we  find  the  expression  for  the  wave  structure  function  appropriate 
to  this  situation.  The  expression  for  the  wave  structure  function  spectrum 
is  obtained  from  the  phase  structure  function  spectrum  by  putting  the 
cosine  term  to  unity  in  Eq.  (B2a); 
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• WW  ifc>m.frd 


(B7)  Fw(<)  - 2ir2  *no(K,0)  ]L  C2(n)  dn  . 

;o 

Performing  the  longitudinal  integration  gives 


(B8) 


Fw(iO  = .033irk2(3Ho)  C2(Hq)  <0^)  + k2 


\ -11/6 
> 


In  the  inertial  subrange,  k2  » (1.077/Lo)2,  the  associated  structure 
function  is 


(B9a)  Dw(p)  = 2 x .033TrkZ (3Ho)C^(Hq)  x 4,1  (1-Jo(kP))k’8/3  d< 


(B9b) 

(BIO) 

Combining 

(Bll) 


, 4 , .033.2  x (6/5)r(V6)  fcZ(3H  )c2(„  )(5/3 
r(ll/6)  25/3  0 n 0 

= 6.88(p/rQ)5/3  . 

Eqs . (B9b)  and  (BIO)  we  have 
rQ  = (6.88/2.91k2(3HQ)C3(Ho))3/5  . 


It  is  interesting  to  note  that  the  value  for  rp  thus  obtained  corresponds 
to  the  plane  wave  value  for  horizontal  propagation  at  the  height  of  the 
receiver  and  over  a path  length  of  only  three  times  the  receiver  height. 


Substituting  for  k2(3H0)C2(Hp)  in  Eq.  (B8)  gives  the  desired 
normalized  form  of  the  phase  structure  spectrum. 


(B12)  F$(k)  = .245  r0‘5/3 


(1  .077/Lo)2  + <2! 


-11/6 
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For  many  purposes  the  receiver  is  much  smaller  than  the  outer  scale, 
(the  outer  scale  being  of  the  same  order  of  magnitude  as  the  height). 
For  such  cases  the  outer  scale  can  be  neglected,  giving 

(B13)  F (k)  = .24b  r '5/3  *"11/3 

5 C 


Equations  (12)  and  (13)  are  the  main  results  of  this  appendix. 
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